AN EFFICIENT MULTIGRID CALCULATION OF THE FAR FIELD MAP 

FOR HELMHOLTZ EQUATIONS 
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Abstract. In this paper we present a new highly efficient calculation method for the far field amplitude patterns 
that arise from scattering problems governed by the oi-dimensional Helmholtz equation. The method is based upon a 
reformulation of the standard real-valued Green's function integral expression for the far field amplitude on a complex 
contour. On this complex contour the scattered wave can be calculated very efficiently using the iterative multigrid 
method, resulting in a fast and scalable calculation of the far field mapping. The full complex contour approach is 
successfully validated on model problems in two and three spatial dimensions. 

Key words. Far field map, Helmholtz problem, multigrid method, complex contour. 

1. Introduction. Scattering problems are of key importance in many areas of science and 
engineering since they carry information about an object of interest over large distances, remote 
from the given target. Consequently, ever since their original statement a variety of applications of 
scattering problems have arisen in many different scientific subdomains. In chemistry and quantum 
physics, for example, virtually all knowledge about the inner workings of a molecule has been 
obtained through scattering experiments 31 . Similarly, in many real-life electromagnetic or acoustic 
scattering problems information about a far away object is obtained through radar or sonar |15j . 
intrinsically requiring the solution of 2D or 3D wave equation. 

Preconditioned Krylov subspace methods are currently among the most efficient numerical al- 
gorithms for the solution of high-dimensional equations, as they exploit the sparsity structure of the 
discretized system of equations and allow for reasonably good scaling with respect to the number of 
unknowns. Indeed, preconditioned Krylov subspace methods are able to solve some symmetric posi- 
tive definite systems in only 0(n) iterations, where n is the number of unknowns in the system |44j . 
However, scattering problem are often described by Helmholtz equations, which after discretization 
lead to highly indefinite linear systems that are notoriously hard to solve using the current genera- 
tion of iterative methods. Moreover, the highly efficient iterative multigrid method [9l fTTT [12l [42l [43] 
is known to break down when applied to these type of problems |17[ 123] . 

Over the past decade significant research has been performed on the construction of good pre- 
conditioners for Helmholtz problems. Recent work includes the wave-ray approach |10j . the idea of 
separation of variables [33], algebraic multilevel methods [8], multigrid deflation [40] and a trans- 
formation of the Helmholtz equation into an advection-diffusion-reaction problem [24]. In 2004 the 
Complex Shifted Laplacian ( CSL ) preconditioner was proposed by Erlangga, Vuik and Oosterlee 
[2"ffl r2~il |2"2] as an effective Krylov subspace method preconditioner for Helmholtz problems. The 
key idea behind this preconditioner is to formulate a perturbed Helmholtz problem that includes a 
complex valued wave number. Given a sufficiently large complex shift, this implies a damping in 
the problem, thus making it solvable using multigrid in contrast to the Helmholtz problem with real 
valued wave numbers. By introducing the complex shifted problem as a preconditioner, the resulting 
Krylov method has advantageous spectral properties, leading to a reasonable convergence rate. The 
concept of CSL has been further generalized in a variety of papers among which [21 IT81 [T^l 132] , 

Recently a variation on the Complex Shifted Laplacian scheme by the name of Complex Stretched 
Grid (CSC) was proposed in [361 137] . introducing a complex valued grid distance instead of a 
complex valued wave number in the preconditioning system. It was furthermore shown in [35] that 
the resulting Krylov subspace method has very similar convergence properties. Indeed, the CSC 
preconditioner can be shown to be generally equivalent to the CSL scheme, and can thus be solved 
equally efficiently using multigrid. 

The choice of a sufficiently large complex shift parameter, denoted in the literature by /3, is vital 
to the stability of the multigrid solution method. The general rule of the thumb for the choice of 
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the complex shift suggested in the literature is to take j3 = 0.5 [3TJ[35]. This experimental guideline 
was recently confirmed through a rigorous LFA analysis in |f 6j . proving multigrid to be generally 
stable for shifts j3 larger than 0.5. 

However, despite its overall qualitative performance the CSL/CSG preconditioned Krylov sub- 
space solution method suffers from a significant wave number dependency of the convergence rate 
[35] - Additionally, the convergence rate quickly deteriorates in the presence of evanescent waves in 
the Helmholtz equation. 

This paper focuses on calculating the far field mapping resulting from a Helmholtz scattering 
problem. Typically, the calculation of the far field map is a two step process. First a Helmholtz 
problem with absorbing boundary conditions is solved on a finite numerical box covering the object 
of interest. In the second step a volume integral calculates the angular dependency of the far field 
amplitude with an integral over the Green's function and the numerical solution. This strategy 
was successfully applied to calculate impact ionization in hydrogen [38] and double photo-ionization 
in molecules [46] I47j described by the Schrddinger equation, which in this case translates into a 
6-dimensional Helmholtz problem. 

The absorbing boundary conditions used in this paper are based on the principle of Exterior 
Complex Scaling (ECS) that was introduced in the I970's [TJ [H HI] and is frequently used in appli- 
cations. This method is equivalent to a complex stretching implementation of PML [14] . 

In this paper we propose a new method for the calculation of the far field map. The method 
reformulates the integral over the Green's function on a complex contour. This modified approach 
requires the solution on the Helmholtz equation on a complex contour. It is shown that the latter 
problem is equivalent to a Complex Shifted Laplacian problem that can be solved very efficiently by 
using a multigrid method. To validate our approach, the method is successfully illustrated on both 
2D and 3D Helmholtz and Schrddinger equations for a variety of discretization levels. 

The outline of the article is the following. In Section 2 we briefly define the notation and 
terminology used throughout the text, and we demonstrate the standard calculation of the far field 
map for Helmholtz type scattering problems. In the second part of this section we introduce an 
alternative way of calculating the far field mapping based upon a reformulation of the integral 
over a complex contour, for which the corresponding Helmholtz system is very efficiently solved 
iteratively. The new technique is validated on a variety of model scattering problems in both 2D 
and 3D in Section 3, where it is found to yield a very fast and scalable far field map calculation 
method. Finally, conclusions are drawn in Section 4. 

2. The Helmholtz equation and the far field map. In this section we introduce the 
general notation used throughout this text and we illustrate the derivation of the far field solution 
and calculation of its amplitude from a Helmholtz-type scattering problem. 

2.1. Derivation of the far field mapping. The Helmholtz equation is a simple mathematical 
representation of the physics behind a wave scattering at an object defined on a compact support 
area O located within a domain f2 C R d . The equation is given by 



with dimension d > 1, where A is the Laplace operator, / designates the right hand side or source 
term, and k is the (spatially dependent) wave number, representing the material properties inside 
the object of interest. Indeed, the wave number function k is defined as 



where k £ K is a scalar constant denoting the wave number outside the object of interest. The 
scattered wave solution is given by the unknown function u. Throughout the text we will use the 
following convenient notation 



(-A- k 2 {x))u{x) =f(x) on fief 
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such that k 2 (x) = kg (1 + x( x ))- Note that the function \ is trivially zero outside the object of 
interest O where the space-dependent wave number k 2 (x) is reduced to ko. Defining the incoming 
wave as Ui n {x) — e lkoT1 ' x , where r\ is the unit vector that defines the direction of the incoming wave, 



the right-hand side f(x) is typically given by k 2 x{ x )u in {x) . Reformulating (2.1 1, we obtain 



(-A - k 2 {x)) u{x) = k 2 x(x)u in (x) for x e f2. (2.4) 

This equation is typically formulated on the domain Q with outgoing wave boundary conditions on 
<9f2. The above equation can in principle be solved in a numerical box (i.e. a discretized subset of £1) 
covering the support of Xj with absorbing boundary conditions along all boundaries. Let us assume 



that the numerical solution satisfying (2.4) on this box has been calculated and is denoted by u 



N 



In order to calculate the far field scattered wave pattern the above equation is reorganized as 
(-A - k 2 ) u(x) = k 2 X (x) (u m {x) + u{x)) . (2.5) 
Note that we can replace the function u(x) in the right hand side of this equation with the nu- 



merical solution u (x) obtained from equation (2.4 1. In doing so, the above equation becomes an 



inhomogeneous Helmholtz equation with constant wave number 

(-A - k 2 ) u{x) = g{x) for x e R d , (2.6) 

where the short notation g(x) = koX{ x )( u in{ x ) +u N (x j) is introduced for readability and notational 
convenience. The above equation can easily be solved analytically using Green's function G(x,x'), 
i.e. 

u(x)= I G{x,x')g(x')dx'. (2.7) 



Since the function g is only non-zero inside the numerical box that was used to solve equation (2.4 1, 
the above integral over K d can be replaced by a finite integral over i7 

u(x) = [ G{x,x')klx{x') (win(x') +u N (x')) dx' for x € M. d . (2.8) 
Jn 

This expression allows us to calculate the scattered wave solution u in any point x € M. d outside the 
numerical box, using only the information inside the numerical box. 



Given the integral expression (2.8), the asymptotic form of the Green's function can be used to 
compute the far field mapping of the scattered wave u. In the following this will be illustrated for a 
2D model example where the Green's function is given explicitly by 

G(x,x')= i -H^(k \x-x'\), (2.9) 

where i represents the imaginary unit and Hq is the 0-th order Hankel function of the first kind. 
Note that an analogous derivation can be performed in 3D, where we mention for completeness that 
the Green's function is given by 

ik \x-x'\ 

G(x,x') = --. -. (2.10) 

An\x-x'\ 

To calculate the angular dependence of the far field map, the direction of the unit vector a 
is introduced that is in 2D defined by a single angle a with the positive horizontal axis, i.e. a = 
(cosa,sina) T . Rewriting the spatial coordinates x in polar coordinates as x = (pcosa, psina) T 
the asymptotic form of the Green's function for \x\ ^> 1 (p — > oo) is given by 

l -H^\k Q \x-x'\) = -^\\^™^—^/ kop z~ lk - r '' K ''' 
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where we have used the fact that the Hankel function Hq 1 ^ is asymptotically given by 

H oH r ) = \f^, c *v( i ( r - \)) > reM ' r>1 - ( 2 - 12 ) 

This leads to the following asymptotic form of the 2D scattered wave solution 

u(p,a) = -A^-^^^L [ e- ik ° x ' a g{x')dx', (2.13) 



4 V 7T y/kop , 



n 



for p — > oo. The above expression is called the 2D far field wave pattern of u, with the integral being 
denoted as the far field ( amplitude ) map 

F{a)= f e- lkoX ' a g{x')dx'. (2.14) 
Jn 

Note that the value of the integral only depends on the direction a (or, in 2D, on the angle a) and 



the wave number fco- Expression (2.13) readily extends to the (i-dimensional case, where it holds 
more generally that 

lim u(p, a) = D(p)F(a), a 6 R d , (2.15) 

p— >oo 



for a function D(p) which is known explicitly and a far field mapping F(a) given by (2.14). This 
far field mapping is in fact a Fourier integral of the function g. 

Summarizing, we conclude that the calculation of the far field wave pattern of the scattered 
wave u consists of two main steps. First, one has to solve a Helmholtz equation with a spatially 



dependent wave number on a numerical box with absorbing boundary conditions as in (2.4). Once 
the numerical solution is obtained, it is followed by the calculation of a Fourier integral (2.14) 
over the aforementioned numerical domain. The main computational bottleneck of this calculation 
generally lies within the first step, since one requires a suitable (iterative) method for the solution 
of a high dimensional indefinite Helmholtz system with absorbing boundary conditions. 

Note that the statement of the far field mapping presented in this section relies on the fact that 
the object of interest function \ IS compactly supported. In particular, this is used when computing 



the numerical solution u N to equation (2.4) on a bounded numerical box that covers the support 
of x- The above reasoning can however be readily extended to the more general class of analytical 
object functions \ that vanish at infinity, i.e. \ <E V where V = {/ : R d — > K analytical | Ve > 
0, 3K C K d compact, Vx € M. d \ K : \f{x)\ < e}. Indeed, due to the existence of smooth bump 
functions [271 130) , functions with compact support can be shown to be dense within the space of 
functions that vanish at infinity. Consequently, every analytical function \ that vanishes at infinity 
can be arbitrarily closely approximated by a series of compactly supported functions {xn}n- This 
in turn implies that the corresponding solutions {u^ } n on a limited computational box can be 
arbitrarily close to the solution of the Helmholtz equation generated with the analytical object of 
interest x G V ■ Intuitively, this means that if x is analytical but sufficiently small everywhere 
outside O, the computational domain may be retricted to a numerical box covering O as though \ 



was compactly supported. Hence, the far field mapping (2.14) is well-defined for analytical functions 



X that vanish at infinity. This observation will prove particularly useful in the next section. 

2.2. Calculation on a complex contour. In this paragraph we will illustrate how the integral 



(2.14) can be reformulated on a complex contour and why this is useful in terms of numerical 
computation. First, we note that the integral can be split into a sum of two contributions F(ot) = 
I\ + 12 with 

h = [ e- lkoX - a x(x)u ln {x)dx and h = \ e- lk ° x - a X {x)u N (x)dx, (2.16) 

Calculation of first integral 1\ is generally easy, since it only requires the expression for the incoming 
wave, which is known analytically. The second integral however requires the solution of the Helmholtz 



4 



■-••••Tf 



M{z} 



Fig. 2.1. Schematic representation of the complex contour for the far field integral calculation illustrated in ID. 
The full line represents the real domain Q, the dotted and dashed lines represent the subareas Z\ = {xe l t : x £ fl} 
and Z2 = {be 1 ® : b £ c9f2, 9 £ [0,7]} of the complex contour respectively. 

equation on the numerical box, which is known to be notoriously hard to obtain using iterative 
methods. In particular, the highly efficient multigrid solution method is unable to solve these type 
of indefinite Helmholtz equations due to instability in both the coarse grid correction and relaxation 
scheme. This divergence is due to close-to-zero eigenvalues of the discretized operator on some 
intermediate multigrid levels |17j . 

However, if both u and \ are analytical functions the integral can be calculated over a complex 
contour rather than the real axis as follows. Let us define a complex contour along the rotated real 
domain Z\ = {z g C | z = xe 11 : x € CI}, where 7 is a fixed rotation angle, followed by the cu rved 
segment Z2 = {z e C | z = be 10 : b £ dfl, < 9 < 7}, as presented schematically on Figure 2.1 
The integral I2 can then be written as 

I 2 = f er lkaZ - a x(z)u N (z)dz+ f e- lkoZ - a x(z)u N (z)dz. (2.17) 

The second term in the above expression however vanishes, as the function \ is P er definition zero 
everywhere outside the object of interest O, thus notably in all points z € Z^. Hence one ultimately 
obtains 

I 2 = I e- ikoz a X (z)u N (z)dz= [ e~' koe "' x - a x(xe i ' 1 )u N (xe l ' r )e i ' 1 dx. (2.18) 

Note that for < 7 < 7r/2 the exponential of xe 11 is increasing in all directions. At the same 
time the scattered wave solution u , which consists of outgoing waves on the complex domain Z\, 
is decaying in all directions. Additionally, the function x is presumed to have a bounded support 
making the above integral calculable on a limited numerical domain. 



Expression (2.181 for the integral I2 indicates that the far field map can (at least partially) be 
computed over the full complex contour Z\, i.e. a rotation of the original real domain fl over an 
angle 7 in all spatial dimensions. The advantage of this approach is that we only need the val ue of 



u N evaluated along this complex contour; thus we now have to solve the Helmholtz equation (2.4) 
on a complex contour. On this contour it is a damped equation which is much easier to solve that 
the Helmholtz equation along the real axis. Indeed, given a sufficiently large value of 7, it has been 
shown in the literature [551 13S] that the multigrid scheme is a very effective solution method for the 
Helmholtz equation on a complex domain. 

2.3. Solving the Helmholtz equation on a complex contour. We now show that the 
Helmholtz problem on the complex domain Z\ is similar to a complex shifted Laplacian system |20j , 
and can thus be solved very efficiently using a multigrid solver. Consider the Helmholtz problem 
with a complex shifted wave number 

(-A- {l + i(3)k 2 (x))u(x) = f(x) (2.19) 

with Dirichlet boundary conditions u(x\gci) = and a complex shift parameter j3 € R. After finite 
difference discretization on a d-dimensional Cartesian grid with fixed grid distance h in every spatial 
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Fig. 2.2. Real grid with ECS boundaries vs. full complex grid. 



dimension, one typically obtains a linear system 

- (^L + {I + ip)k 2 ^j u h =b h (2.20) 

where L is the matrix operator expressing the stencil structure of the Laplacian. In 2D, for example 
L = kron(J, diag(— 1, 2, —1) + kron(diag(— 1, 2, —1), /), where the size of L intrinsically depends on 



h. After dividing both sides in linear system (2.20) by (l+i/3), we immediately obtain the equivalent 
system 

' • L + k2 ) u ^-^~ 5 , (2.21) 



(l + iP)h 2 J 1 



which is identical to the discretization of the original Laplacian with grid distance h = y/T+Tph. 
This scheme is known as Complex Stretched Grid, and it was shown in |36j to yield exactly the same 
Krylov convergence compared to Complex Shifted Laplacian when both are used as a preconditioner 
for a general Krylov method. 



It is known from the literature [THIHO] that problem (2.20), or equivalently (2.21 ), can be solved 
efficiently with multigrid for values of the complex shift f3 > 0.5. Note that this requirement is based 
on a multigrid cycle with standard weighted Jacobi or Gauss-Seidel smoothing. This rule of thumb 
can easily be translated into an angle 7 for the complex scaling. Writing (1 + if3) = pexp(i(p) with 
p = yl + j3 2 and tp = arctan/3, one readily obtains 

h = y/l + i/3 h = x /pexp(iip/2)h (2.22) 

Consequently, as the shift (3 is required to be larger than 0.5, the grid rotation angle 7 = ip/2 must 
satisfy 

7 > arCta 2 n(0 ' 5) = 0.2318 « 13.28° (2.23) 

Note that when substituting the standard multigrid relaxation schemes like w- Jacobi or Gauss-Seidel 
by a more robust iterative scheme like e.g. GMRES(m), the rotation angle 7 may be chosen even 
smaller, up to a minimum of approximately 9.5° (see |34j ) . 

In this paper we have chosen to link the grid rotation angle 7 to the standard ECS absorbing 
layer angle 9, see Figure [2~2| This is in no way imperious for the functionality of the method, but it 
appears quite naturally from the fact that both angles perturb (part of) the grid into the complex 
plane. Suppose the ECS boundary layer measures one quarter of the length of the entire real domain 
in every spatial dimension, which is a common choice, we readily derive that the relation between 
the rotation angle 7 and the ECS angle 9 is given by 

/ sin9 \ . „ 

7 = arctan . (2.24) 



2 + cos 6», 

Table |2.1| shows some standard values of the ECS angle 9 and corresponding 7 values according to 



(2.24). Note that for a multigrid scheme with ui- Jacobi or Gauss-Seidel smoothing to be stable, 9 



should be chosen no smaller than 7r/4, see (2.23). Using the more efficient GMRES(3) method as a 
smoother, the ECS angle can be chosen somewhat smaller, i.e. an angle around 9 = ir/6 suffices to 
guarantee a stable multigrid solution. 



9 (rad.) 


tt/8 


tt/6 


tt/5 


7r/4 


tt/3 


7 (deg.) 


7.5° 


9.9° 


11.8° 


14.6° 


19.1° 



Table 2.1 

iJCS 1 angle 8 and corresponding rotation angle 7 /or the full complex grid. Values based on SK 



3. Numerical results for 2D and 3D Helmholtz problems. In this section, we validate 
the theoretical result presented above by a number of numerical experiments in both two and three 
spatial dimensions. It will be shown that the proposed method results in a very fast and wave 
number-independent solution method for the scattered wave system, hence yielding a remarkably 
efficient method for the calculation of the far field mapping. 



The model problem used throughout this section is a Helmholtz equation of the form (2.4) 
with k 2 (x) = A;q (1 + x( x ))- The equation is discretized on a n d -point uniform mesh covering a 
square numerical domain £1 — [—20, 20] d using second order finite differences. In the 2D case the 
space-dependent wave number is defined as 

X (x, y) = -1/5 (e-(* 2 +(^) 2 ) + e -^+(v+^ , fa y) e [_ 2 0, 20] 2 , (3.1) 

i.e. the object of interest takes the form of two circular point-like objects with mass concentrated 
at the Cartesian coordinates (0,-4) and (0,4) (see Figure 3.1). For the 3D model problem, the 
following straightforward extension of the object is used 

X (x, y,z) = -1/5 ( e -(- 2 +fe- 4 ) 2 + z2 ) + e -(x 2 +(y+D 2 +z 2 )^ ) fa y ^ z) e [_ 20 , 20] 3 , (3.2) 



representing two spherical point-like objects in 3D space (see Figure 3.2). The incoming wave 
scattering at the given object is defined by 

u in (x) = e !fcnI ^, x e n, (3.3) 

where rj is the unit vector in the x-direction. 



Figure 3.1 illustrates the theoretical re sult presented in Section 2. The above 2D Helmholtz 
model problem with wave number given by (3.1) is solved for ti" using respectively a standard LU 



factorization method on the real domain Q with ECS complex boundary layers {6 = 7r/4) along 
the domain boundary 9f2, and a series of multigrid V(l,l)-cycles with w-Jacobi smoothing on the 
full complex domain (7 sa 14.6°) up to residual tolerance le-6. The standard multigrid intergrid 
operators used in this work are bilinear interpolation and full weighting restriction. The mo duli 
of the wave number function \ (top) and the resulting solution u N (mid) are shown on Figure 3.1 



for both methods. Note how the solution u N on the full complex contour is indeed heavily damped 
compared to the solution on the real domain. Consequently, using the numerical solution u N , the 2D 



far field map integral (2.14) can be calculated using any numerical integration scheme over the real 



or complex domain respectively. The resulting far field mapping F(a) is shown as a function of the 



angle a on Figure 3.1 (bottom). One observes that the mapping is indeed identical when calculated 
on the real and complex domain, conform with the theoretical results. However, the computational 
cost of the real-domain method for calculation of the far field map is reduced significantly by the 
ability to apply multigrid to the equivalent complex scaled problem. 

In Table |3.1| convergence results are shown for the solution of the 3D scattered wave equation 



(2.4) using a series of multigrid V(l,l)-cycles on various grid sizes. Note that the multigrid method 
scales perfectly as a function of the number of grid points, as doubling the number of grid points in 
every spatial dimension does not increase the number of V-cycles required to reach a fixed residual 
tolerance of le-6. This is a standard result from multigrid theory. Additionally and more im- 
portantly, remarkable fc-scalability is measured for the multigrid solution method on the complex 
contour. Indeed, the multigrid convergence factor (and thus the corresponding work unit load re- 
quired to solve the problem up to a given tolerance) is almost fully independent of the wave number 
kg, as can be observed from the table. Note that from a physical-numerical point of view it is only 
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Fig. 3.1. Top: 2D object of interest \x( x )\ given by ^3. IV . Mid: solution to the Helmholtz problem \2.J$ (in 
modulus) on a n x X n y = 256 X 256 grid solved using LU factorization (left) on a double ECS contour with 8 = tt/4, 
and using a series of multigrid V-cycles (right) with uj-Jacobi smoother on the corresponding full complex contour 
up to a residual reduction tolerance of le-6. Bottom: resulting 2D Far field maps F(a) calculated following \2. 1J$. 



meaningful to consider discretizations satisfying the k a h < 0.625 criterion for a minimum of 10 grid 
points per wavelength, cfr. [6 , for which the corresponding values are designated in Table 3.1 by a 
bold typesetting. 

Ultimately, the computed scattered wave solution on the complex domain can again be used to 
calculate the far field integral (2.14). The resulting 3D far field mapping for the model problem with 
ko = 1 is plotted in Figure 3.2 The left hand side panel shows an iso-surface visualization of the 
3D object of interest x(x) given by (3.2). On the right panel a spherical projection of the resulting 
3D far field mapping is shown. The color hue indicates the value of the far field amplitude in each 
outgoing direction. 

Note that the calculation of the scattered wave solution can be optimized even further by 
considering the Full Multigrid (FMG) scheme. This is a nested iteration of standard V-cycles, where 
on each level a series of V(l,l)-cycles is used to approximately solve the error equation and supply 
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Fig. 3.2. Left: 3D object of interest \x(x)\ given by \3. 2[ ) . Shown are the = c isosurfaces for c = 0, 

le-2, le-10 and le-100. Right: 3D Far field map, resulting from Helmholtz problem \2.4\ with ko = 1 solved on a 
n x X n y X n z = 64 X 64 X 64 full complex grid with 8 = 7r/6 (j ~ 9.9°) using a series of multigrid V-cycles with 
GMRES(S) smoother up to residual reduction tolerance le-6. 







16 3 


32 3 


64 3 


128 3 


256 3 




1/4 


10 (10) 


9 (59) 


9 (560) 


9 (4456) 


9 (35165) 




0.24 


0.20 


0.21 


0.20 


0.20 




1/2 


12 (12) 


10 (63) 


10 (611) 


10 (4937) 


9 (35405) 




0.31 


0.24 


0.22 


0.23 


0.21 


fc 


1 


7(8) 


13 (83) 


11 (691) 


10 (4899) 


10 (38975) 




0.13 


0.32 


0.27 


0.24 


0.24 




2 


2(4) 


8 (54) 


13 (809) 


11 (5418) 


10 (38051) 






0.01 


0.14 


0.33 


0.27 


0.24 




4 


1 (3) 


2 (17) 


7 (457) 


13 (6337) 


11 (41848) 






0.01 


0.01 


0.12 


0.33 


0.26 



Table 3.1 

3D Helmholtz problem \2.J$ solved on a full complex grid with 6 = 7r/6 using a series of multigrid V(l,l)-cycles 
with GMRES(3) smoother up to residual reduction tolerance le-6. Displayed are the number of V-cycle iterations, 
number of work units and average convergence factor for various wave numbers ko and different discretizations. 1 
WU is calibrated as the cost of 1 V(l,l)-cycle on the 16 3 -points grid ko = 1/4 problem. Discretizations respecting 
the kgh < 0.625 criterion for a minimum of 10 grid points per wavelength are indicated by a bold typesetting. 



a corrected initial guess for a finer level by interpolating the corresponding coarse grid solution. 



Table 3.2 shows convergence results for the solution of the 3D scattered wave equation (2.4) 



using an FMG scheme. The setting is comparable to that of Table as a residual reduction 
tolerance of 10 -6 is imposed for each wavenumber and at every level of the FMG cycle, yielding a 
fine n x x n y x n z = 256 3 grid residual of order of magnitude 10~ 9 . Note that the number of V-cycles 
performed on each level in the FMG cycle is decaying in function of the growing grid size due to 
the increasingly accurate initial guess, resulting in a relatively small number of V-cycles (5-9) to be 
performed on the finest level. Consequently, the number of work units (and thus the computational 
time) required to reach the designated residual reduction tolerance is significantly lower than the 
work unit load of the pure V-cycle scheme displayed in Table |3.1| 

Timing and residual results from a standard FMG sweep performing only one V(l,l)-cycle 
on each level on the 3D Helmholtz scattering problem with a moderate wave number ko = 1 are 
shown in Table 3.3 for different discretizations. Note that timings were generated using a basic 
non-parallelised Matlab code, using only a single thread on a simple midrange personal computer 
and taking less than 8 minutes to solve a 3D Helmholtz problem with 256 gridpoints in every spatial 
dimension. 
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Tlj^ X Tly X Tbg 


16 3 


32 3 


64 3 


128 3 


256 3 




1/4 


8 (11) 


O ^OZ ) 


o (^004 J 










1 77p_Q 


9 4fip-Q 


2 RSp-Q 






1/2 


9 (12) 


8 (68) 


6 (452) 


6 (3392) 


6 (30215) 




1.27e-8 


1.87e-9 


3.37e-9 


1.30e-9 


1.25e-9 




1 


5(8) 


9 (68) 


8 (572) 


7 (4013) 


6 (30747) 




1.33e-8 


1.51e-8 


4.07e-9 


1.76e-9 


3.66e-9 




2 


1(5) 


5 (43) 


9 (600) 


8 (4456) 


7 (36367) 






5.99e-13 


1.18e-8 


1.91e-8 


5.28e-9 


2.67e-9 




4 


1(4) 


1(18) 


5 (357) 


9 (5038) 


8 (39038) 






8.90e-20 


2.86e-13 


5.19e-9 


1.97e-8 


4.65e-9 



Table 3.2 

3D Helmholtz problem \2.4\ solved on a full complex grid with 9 = n/6 using an FMG cycle with GMRES(3) 
smoother up to residual reduction tolerance of le-6. Displayed are the number of V-cycle iterations on the designated 
finest grid, number of work units and resulting residual norm for various wave numbers ko and different discretiza- 
tions. 1 WU is calibrated as the cost of 1 V(l,l)-cycle on the 16 3 -points grid kg = 1/4 problem. Discretizations 
respecting the koh < 0.625 criterion for a minimum of 10 grid points per wavelength indicated by a bold typesetting. 



X fly X Tt 2, 


16 3 32 3 64 3 128 3 256 3 


CPU time 
IM| 2 


0.20 s. 0.78 s. 6.24 s. 53.3 s. 462 s. 
3.3e-5 7.9e-5 2.7e-5 l.le-5 4.6e-6 



Table 3.3 

3D Helmholtz problem with wave number k(, = 1 solved on a full complex grid with 9 = tt/6 using one 

FMG-cycle with GMRES(3) smoother. Displayed are the CPU time (in s.) and the resulting residual norm for 
various discretizations. System specifications: Intel® Core™ i7-2720QM 2.20GHz CPU, 6MB Cache, 8GB RAM. 



4. Conclusions and discussion. In this paper wc have developed a novel highly efficient 
method for the calculation of the far field map resulting from c£-dimensional Helmholtz scattering 
problems where the wave number is an analytical function. Our approach is based on the refor- 
mulation of the classically real-valued Green's function volume integral for the far field map to an 
equivalent volume integral over a complex valued domain. 

The advantage of the proposed reformulation lies in the scattered wave solution of the Helmholtz 
problem on a complex domain, which for high dimensional problems can be calculated efficiently 
using a multigrid method. Indeed, the reformulation of the Helmholtz forward problem on the full 
complex contour is shown to be equivalent to a Complex Shifted Laplacian problem, where multigrid 
has been proven in the literature to be a fast and scalable solver. However, whereas the Complex 
Shifted Laplacian was previously only used as a preconditioner for highly indefinite Helmholtz prob- 
lems, the complex-valued far field map calculation proposed within this paper effectively allows for 
multigrid to be used as a solver. 

The functionality of the method is validated on 2D and 3D model Helmholtz problems. It is 
confirmed that the values of the far field mapping calculated on the full complex grid exactly matches 
the values of the classical real-valued integral. Furthermore, the number of multigrid iterations is 
shown to be largely wave number independent, yielding a fast overall far field map calculation. 

One area of scientific computing where the proposed technique might be particularly valuable 
is in the numerical solution of quantum mechanical scattering problems. These are generally high- 
dimensional scattering problems where the wave number is indeed an analytical function and where 
6D or 9D problems are common. 

Finally we note that a number of modifications can be made to improve the efficiency of the 
method even further, e.g. choosing the complex contour for the integral based on a steepest descent 
scheme as proposed in [2T>] . 
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